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Constructing a bivariate distribution function with given marginals 
and correlation: application to the galaxy luminosity function 
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ABSTRACT 

We show an analytic method to construct a bivariate distribution function (DF) with given 
marginal distributions and correlation coefficient. We introduce a convenient mathematical 
tool, called a copula, to connect two DFs with any prescribed dependence structure. If the 
correlation of two variables is weak (Pearson's correlation coefficient \p\ < 1/3), the Farlie- 
Gumbel-Morgenstern (FGM) copula provides an intuitive and natural way for constructing 
such a bivariate DF. When the linear correlation is stronger, the FGM copula cannot work any- 
more. In this case, we propose to use a Gaussian copula, which connects two given marginals 
and directly related to the linear correlation coefficient between two variables. Using the cop- 
ulas, we constructed the BLFs and discuss its statistical properties. Especially, we focused on 
the FUV-FIR BLF, since these two luminosities are related to the star formation (SF) activ- 
ity. Though both the FUV and FIR are related to the SF activity, the univariate LFs have a 
very different functional form: former is well described by the Schechter function whilst the 
latter has a much more extended power-law like luminous end. We constructed the FUV-FIR 
BLFs by the FGM and Gaussian copulas with different strength of correlation, and exam- 
ined their statistical properties. Then, we discuss some further possible applications of the 
BLF: the problem of a multiband flux-limited sample selection, the construction of the SF 
rate (SFR) function, and the construction of the stellar mass of galaxies (M*)— specific SFR 
(SFR/M*) relation. The copulas turned out to be a very useful tool to investigate all these 
issues, especially for including the complicated selection effects. 

Key words: dust, extinction - galaxies: evolution - galaxies: luminosity function, mass func- 
tion - infrared: galaxies - method: statistical - ultraviolet: galaxies 
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1 INTRODUCTION 



A luminosity function (LF) of galaxies is one of the f undamental tools to describe an d expl ore the distribution of luminous matter in 
the Universe, (see, e.g. iBinggeli. Sandaee. & Tammannl 1 19881 : Lin et alJ Il99fj : iTakeuchil l200d iTakeuchi et al. 200(j; iBlanton et ail 12001 



de Lapparent et al]2003l : lwillmer et al.l2006h . Up to now, studies on the LFs have been rather restricted to a univariate one, i.e. LFs based on 



a single selection wavelength band. However, such a situation is drastically changing in the era of large and/or deep Legacy surveys. Indeed, 
a vast number of recent studies are multiband-oriented: they require data from various wavelengths from the ultraviolet (UV) to the infrared 
(IR) and radio bands. A bivariate LF (BLF) would be a very convenient tool in such studies. However, to date, it is often defined and used in a 
confused manner, without careful consideration of complicated selection effects in both bands. This confusion might be partially because of 
the intrinsically complicated nature of multiband surveys, but also because of the lack of proper recipes to describe a BLF. Then, the situation 
will be remedied if we have a proper analytic BLF model. 

However, it is not a trivial task to determine the corresponding bivariate function from its marginal distributions, if the distribution is 
not multivariate Gaussian. In fact, there exist infinitely many distributions with the same marginals because the correlation structure is not 
specified. In general astro nomical applications (not only BLFs), for i nstance, a biv ariate distribution is often obtained by either an ad hoc 
or a heuristic manner (e.g. Choloniew skl 1985 : Chapman et al. 20031 : Schafei 2007 1. though these methods are quite well designed in their 
purposes. Further, analytic bivariate distribution models are often required to interpret the distributions obtained by nonparametric methods 
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Driver ct al 



2006)]. For such purposes, a general method to construct a bivariate distribution 



(e.g. ICross & Driveilbood ; iBall et aljliooj ; 
function with pre-defined marginal distributions and correlation coefficient is desired. 

In econometrics and mathematical finance, such a function has been commonly used to analyze two covariate random variables. This 
is called "cop ula". Especially in a bivar iate context, copulas are useful to define nonparametric measures of dependence for pairs of random 
variables (e.g. lTrivedi & Zimmenl2005l) . In astrophysics, however, it is o nly recently that copulas attract researchers' attention and are not 
very widely known yet (still only a handful of astrophysical applications: iBenabed et al 1 120091 : Jiang et al ] |2009l : lKoeiJ[2009l ; I Scherrer et al 



2010). Hence the usefulness and limitations of copulas are still not well understood in the astrophysical community. 



In this paper, we first introduce a relatively rigorous definition of a copula. Then, we choose two specific copulas, the Farlie-Gumbel- 
Morgenstern (FGM) copula and the Gaussian copula to adopt for the construction of a model BLF. Both of them have an ideal property that 
they are explicitly related to the linear correlation coefficient. Though, as we show in the following, the linear correlation coefficient is not a 
perfect measure of the dependence of two quantities, this is the most familiar and thus fundamental statistical tool for physical scientists. We 
focus on the far-infrared (FIR)-far-ultraviolet (FUV) BLF as a concrete example, and discuss its properties and some applications. 

This paper is organized as follows: in Section[2]we define a copula and present its dependence measures. We also introduce two concrete 
functional forms, the Farlie-Gumbel-Morgenstern copula and the Gaussian copulas. In Section|3] we make use of these copulas to construct a 
BLF of galaxies. Especially we emphasize the FIR-FUV BLF. We discuss some implications and further applications in Section|4] Section|5] 
is devoted to summary and conclusions. In Appendix [A] we show an iterated extension of the FGM copula. We present statistical estimators 
of the dependence measures of two variables in Appendix|B]to complete the discussion. 

Throughout this paper, we adopt a cosmological model (h, Qmo, Hao) = (0.7,0.3,0.7) (h = //o/100[km s" 1 ]) unless otherwise 
stated. 



2 FORMULATION 
2.1 Copula 

As we discussed in Introduction, there is an infinite degree of freedom to choose a dependence structure of two variables with given marginal 
distribution. However, very often we need a systematic procedure to construct a bivariate distribution function (DF) of two variable^} Copulas 
have a very desirable property from this point of view. In short, copulas are functions that join multivariate DFs to their one-dimensional 
marginal DFs. However, this statement does not serve as a definition. We first introduce its abstract framework, and move on to a more 
concrete form which is suitable for the aim of this work (and of many other physical studies). 
Before defining the copula, we prepare some mathematical concepts in the following. 

Definition 1. Let Si and S 2 be nonempty subsets of R [a union of real number R and (—00,00)]. Let H be a real function with two 
arguments (referred to as bivariate or 2-place) such that Dom H = Si X S 2 . Let B — [xi, x 2 ] x [yi , y 2 ] be a rectangle all of whose vertices 
are in Dom H . Then, the H -volume of B is defined by 

V H (B) = H(x 2 ,y 2 ) - H(x 2 ,yi) - H(x 1 ,y 2 ) + H(xi, yi ) . (1) 

Definition 2. A bivariate real function H is 2-increasing if Vh ^ for all rectangles B whose vertices lie in Dom H. 

Definition 3. Suppose Si has a least element ai and S 2 has a least element a 2 . then, a function H: Si x S 2 — > R is grounded if H(x, a 2 ) — 
and H(ai,y) — for all (x, y) in Si x S 2 . 

With these concepts, we now define a two-dimensional copula. 

Definition 4. A two-dimensional copula (or shortly 2-copula) is a function C with the following properties: 

(i) DomC* = [0,1] x [0,1]; 

(ii) C is grounded and 2-increasing; 

(iii) For every u and v in [0, 1], 

C(u, 1) = u and C(l, v) = v . (2) 
It may be useful to show there are upper and lower limits for the ranges of copulas, which is given by the following theorem. 
Theorem 1. Let C be a copula. Then, for every (u, v) in Dom C, 

max(it + v — 1, 0) ^ C(u, v) ^ min(M, v) . (3) 

Often the notations W(u,v) = max(u + v — 1,0) and M(u,v) = min(u,v) are used. The former and the latter are referred to as 
Frechet-Hoeffding lower bound and Frechet-Hoeffding upper bound, respectively. 

1 In this work, we use a term DF for a cumulative distribution function (CDF). We use a term probability density function (PDF) to avoid confusion with the 
term used in physics "distribution function" which stands for a Radon-Nikodym derivative of a DF. 
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Figure 1. The Frechet-Hoeffding lower and upper bounds. Left and right panels show the lower and upper bounds, respectively. The top panels represent the 
contours of W(u, v)) and M(u, v). 



The Frechet-Hoeffding lower and upper bounds are illustrated in Fig.[T] Any copula has its value between W(u, v) and M(u, v). 

Any bivariate function which satisfies the above conditions can be a copula. Then, there are infinite degrees of freedom for a set of 
copulas. Using a copula C, we can construct a bivariate DF G with two margins Fi and F 2 as 



G(xi,x 2 ) = C[F(x t ),F 2 (x 2 )] 



(4) 



Howe ver, one may have a natural question: can any bivariate DF be written in the above form? This is guaranteed by Sklar's theorem ( Sklar 
1959b . 



Theorem 2. (Sklar's theorem) Let G be a joint distribution function with margins Ft and F 2 . Then, there exists a copula C such that for all 

Xi, x 2 in R, 

G(x 1 ,x 2 ) = C[F 1 (x 1 ),F 2 (x 2 )] . (5) 
If Fi and F 2 are continuous, then C is unique: otherwise, C is uniquely determined on Range Ft x Range F 2 . 

A comprehensive proof Sklar's theorem is found in e.g., Nelsen (2006). This theorem gives a basis that any bivariate DF with given margins 



is expressed with a form of Equation QJ. Then, finally, the somewhat abstract definition of a copula turns out to be really useful for our aim, 
i.e., to construct a bivariate DF when its marginals are known in some way. 

Up to now, we discuss only bivariate DFs and their copulas. It is straightforward to introduce multivariate DFs as a natural extension of 
the formulation presented here. 



2.2 Copulas and dependence measures between two variables 

The most important statistical aspect of bivariate DFs is their dependence properties between variables. Since the dependence can never given 
by the marginals of a DF, this is the most nontrivial information which a bivariate DF provides. Since any bivariate DFs are described by 
Equation l|5}, all the information on the dependence is carried by their copulas. 

The most familiar measure of dependence among physical scientists (and others) may be the correlation coefficients, especially Pearson's 
product-moment correlation coefficient p. The bivariate PDF of Xj and x 2 , g(xi,x 2 ), is written as 

g{xuX2) = &C[Fi{xi),ft{x 2 )} fl(x 1 )f 2 (X2)=c[F 1 (x 1 ),F 2 (X2)]f 1 (x 1 )f2(x 2 ) (6) 

where fi(xi) and f 2 (x 2 ) are PDF of Ft (a;i) and F 2 (x 2 ), respectively. Then the correlation coefficient is expressed as 

^_ J{xi - xi)(x 2 - x 2 )g(xi,x 2 )dxidx 2 _ j(xi - x~!)(x 2 - x 2 )c[(F!(xi), F 2 (x 2 )]f 1 (xi)f 2 (x 2 )dx 1 dx 2 ^ 

J f(x! - xi) 2 fi(x!)dxi J(x 2 - x 2 ) 2 f 2 (x 2 )dx 2 J /(si - si) 2 /i(xi)da;i J(x 2 - x 2 ) 2 f 2 (x 2 )dx 2 

We should note that p can measure only a linear dependence of two variables. However, in general the dependence of two variables would 

© 0000 RAS. MNRAS 000, 000-000 



4 T. T. Takeuchi 



be not linear at all, and it cannot be a sufficient measure of dependence. Further, more fundamentally, Equation depends not only on the 
dependence of two variables (copula part) but also its marginals fi(xi), f 2(^x2), i- e -, the linear correlation coefficient p does not measure the 
dependence purely. In such a situation, a more flexible and genuine measure of dependence, e.g., Spearman's ps or Kendall's r would be 
more appropriate. Spearman's rank correlation is a nonparametric version of Pearson's correlation using ranked data. The population version 
of Spearman's ps is expressed by copula as 



ps = 12 



/ / uiU2dC'(ui, U2) — 3 
Jo Jo 



= 12/ f C(ui,w 2 )duid«2-3. (8) 
Jo Jo 

The definition of Kendall's tau is more complicated. We define a concept of concordance as follows: when we have pairs of data (xii,x 2 i) 
and (x\j,X2j), they are said to be concordant if xu > xij and X2i > X2j or xu < x\j and X2i < X2j (i.e., (xu — X\j)(x2i — X2j) > 0), 
and otherwise discordant. Let {xu, %2i}i=i,...,n denote a random sample of n observations. There are n C2 pairs (xu, X2i) and (xij , X2j) 
of observations in the sample, and each pair is concordant or discordant. Let n c denote the number of concordant pairs, and rid the number 
of discordant ones. Then, Kendall's r for the sample, t, is defined as 

t = . (9) 

n c + n d 

The population version of r is also expressed in a simple form in terms of copula as 

C(iti, w 2 )dG(«i, U2) — 1 

= 4 / / C(ui, U2)c(ui, «2)d«id«2 — 1 ■ (10) 
Jo Jo 

Note that both Equations {8]l and dlOt are independent of the distributions Fi , F2 and G, but they only show the dependence structure 
described by the copula, unlike the linear correlation coefficient. This is a direct consequence of the non-parametric (i.e., distribution-free) 
nature of these estimators. These are the reasons why both dependence measures are almost always used in the context of copulas in the 
literature. Estimators of ps and r from a sample are found in Appendix 151 



2.3 Farlie-Gumbel-Morgenstern (FGM) copula 



As seen in the discussion above, usefulness of the linear correlation coefficient is quite limited, and distribution-free measures of dependence 
are more appropriate for general joint DFs with non-Gaussian marginals. However, even if it is true, since physicists may cling to the most 
familiar linear correlation coefficient p, a copula which has an explicit dependence on p would be convenient. We introduce two special types 
of copulas with this ideal property. 

For cases where the correlation between two variables is weak, a systematic method ha s been propos ed by iMorgensternl 1 19560 and 
Gumbel 1 1960!) for specific functional forms, and later generalized to arbitrary functions by Farlid Jl960l) . This is known as the Farlie- 
Gumbel-Morgenstern (FGM) distributions after the inventors' names. Although the study of the FGM distributions does not seem tightly 
connected to that of copulas, we will see later that they can be expressed i n terms of the so-called FGM co pula. 

The correlation structure of the FGM distributions was studied bv ISchucanv. Parr. & Boverl 1 119781) . Let Fi(x\) and F 2 (a;2) be the 
(cumulative) distribution functions (DFs) of a stochastic variables x and y, respectively, and let and f 2 (x 2 ) be their probability 

density functions (PDFs). The bivariate FGM system of distributions G(xi, X2) is written as 



G(xi, X2) = F 1 (x 1 )F 2 (x 2 ) {1 + k [1 - Fi(xi)] [1 - F 2 (:T2)]} 



(11) 



where G(a;i, 2:2) is the joint DF of Xi and a 2 . Here k is a parameter rel ated to the correl ation (see below), and to make G(x\, X2) have an 
appropriate property as a bivariate DF, | k\ ^ lis required (for a proof, see ICambanisll977l) . Its PDF can be obtained by a direct differentiation 

of G(xi,X2) as 



g{xi,x 2 ) 



d 2 G 



i:) 2 



' dxxdx 2 X1]X2 
-■fi(xi)Mx 2 ){l + K[2F 1 (xi) 



{Fi(xi)F 2 (x 2 ) [1 + k (1 - Fi(xi)) (1 - F2O2))]} 



dxidx2 

1] [2F 2 (x 2 ) - 1]} . 

From Equation U2t . it is straightforward to obtain its covariance function Cov(a;i, X2) as 
Cov(xi, x 2 ) = k J xj [2Fi(a5i) — 1] /i(a;i)dxi J x 2 [2F 2 (x 2 ) - 1] f 2 (x 2 )dx 2 . 

Then we have a correlation function of two stochastic variables Xi and x 2 , p(xi,x 2 ) [Equation Q] as follows 



p(x 1 ,x 2 ) = 



Co\i(xi, x 2 ) 



C102 



<7l<72 



Xl 



[2Fi(a;i) - 1] f 1 (x 1 )dx 1 J x 2 [2F 2 (x 2 ) - 1] f 2 (x 2 )dx 2 



(12) 



(13) 



(14) 
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where (Ji and 02 are the standard deviations of X\ and X2 with respect to /1 (xi) and f% (x?). It is straightforwardly confirmed that g(xi, X2) 
really has the marginals fi(xi) and f2(x2), by a direct integration with respect to xi or X2- It is also clear that if we want a bivariate PDF 
with a prescribed correlation coefficient p, we can determine the parameter k from Equations ((7} and d!3t . 
Here, consider the case that that both fi(x\) and 72(2:2) are the Gamma distributions, i.e., 



x - a ~ 1 e~ x j/ b 
fj{X]) = 1 b"T(a) 

where r(a;) is the gamma function (J = 1, 2). In this case, after some algebra, p can be written analytically as 



p(x!,x 2 ) = -== , . 

Vab L B(a,a) 

where 



2~2(a-l) 



2 -2(6-l) ■ 



B(b,b) 



B(a,b) 



r(q)r(ft) 

T(a + 6) 



(15) 



(16) 



(17) 



D'Estc 



19811) . Especially when 6 = 1, this corresponds to a bivariate extension of the Schechter function 



This result was obtained by 

i Schechteilll976t) , and we expect some astrophysical applications. 

As we mentioned, the correlation of the FGM distributions is restricted to be weak: indeed, the correlation coefficient cannot exceed 
1/3. Here we prove this. For all (absolutely continuous) F(x), 



x [2F(x) - 1] /(i)dx| = jy (x-x) [2F(x) - 1] /(x)dx| 

^ j (x - xf f(x)dx I [2F(x) - l] 2 f(x)dx = ^ , 



(18) 



where x is the average of x. The second line follows from the Cauchy-Schwarz inequality. From Equations J 1 3b and ifTJ, and the condition 
|«| ^ 1, we obtain \p\ ^ 1/3. 

The copula of the FGM family of distributions is expressed as 



C (Ul, Ul\ K) = UlW2 + KUlU2(l — Wi)(l — U2) 

with — 1 < k < 1. The differential form of the FGM copula follows from Equations ((6) and l |19t 



(19) 



(ui, W2; k) = 1 + — 2ui)(l — 2^2) 



(20) 



2.4 Gaussian copula 



As seen in Section 12.31 though the FGM distribution has one of the most "natural" example of bivariate DF which has an explicit /in- 
dependence, the limitation of the correlation coefficient of the FGM family hampers a flexible application of this DF, though there have been 
many attempts to extend its range of application (see Appendix[A}- Then, the second natural candidate may be a copula related to a bivariate 
Gaussian DF. The Gaussian copula has also an explicit dependence on a linear correlation coefficient by its construction. 
Let 



exp 



*i = / ^(x')dx' 

J —OO 

i>i(x\,X2; p) 



vwa-p 2 ) 



exp 



Xi + x\ — 2pX\X2 



2(1 -P 2 ) 



and 



/ ip i (x' 1 ,x'2)dx' 1 dx' 2 

-OO J —OO 

By using the covariance matrix S 

1 p N 



Equation d23t is simplified as 



exp 



x T S" 1 x 



(21) 
(22) 

(23) 
(24) 

(25) 

(26) 
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10* 10* 10'° 10" 10* 10* 10'° 10 ,a 

60/im luminosity L m [L ] 60/im luminosity Z, 60 [L G ] 



Figure 2. The bivariate luminosity functions (BLFs) constructed with the Farlie-Gumbel-Morgenstern copula, with model luminosity functions (LFs) of 
ultraviolet (UV) and infrared (IR)-selected galaxies. The BLFs are normalized so that integrating it over the whole ranges of L\ and L 2 gives one. From 
top-left to bottom-right, the linear correlation coefficient p = 0.0, 0.1, 0.2, and 0.3, corresponding to n = 0.0, 0.33, 0.67 and 1.0, respectively. The contours 
are logarithmic with an interval A log 0^ 2 ) = 0.5 drawn from the peak probability. 



where x = (xi, X2) T and superscript T stands for the transpose of a matrix or vector. 
Then, we define a Gaussian copula C G (ui, m 2 ; p) as 

The density of C G , c G , is obtained as 

d 2 c G {u u u 2 - P ) _ a 2 * 2 [*r 1 («i),*r l W;p] 



(27) 



C (Ul ,U2\p)- 



duidu2 
ip2{x 1 ,x 2 ; p) 

1 

x/(27r) 2 detS 



du±du2 



exp ( — — x T S x x 



2 . exp ry 



2, exp ri 



1 



VdetS 
1 

VdetS 
1 



exp 
exp 
exp 



x J E _1 x-x J Ix 



x T fS^-flx 



It 



VdetS 

where \E' 1 = [^^(iti), ^^ 1 (u 2 )] T and I stands for the identity matrix. The second line follows from Equation l[6j, 



(28) 
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3 APPLICATION TO CONSTRUCT THE BIVARIATE LUMINOSITY FUNCTION (BLF) OF GALAXIES 
3.1 Construction of the BLF 

We define the luminosity at a certain wavelength band by L = vL v (v is the corresponding frequency). Then the luminosity function is 
defined as a number density of galaxies whose luminosity lies between a logarithmic interval [log L, log L + d log L] : 

4> W (L) = -7T— -y j (29) 
d log L 

where we denote log x = log 10 x and In x = log e x. For mathematical simplicity, we define the LF as being normalized, i.e., 

J <^ (1) (L)dlogL = 1 . (30) 

Hence, this corresponds to a probability density function (PDF), a commonly used terminology in the field of mathematical statistics. We 
also define the cumulative LF as 

«W(L)= f l ° eL (1) (L')dlogL', (31) 

J log Zj m i n 

where L mln is the minimum luminosity of galaxies considered. This corresponds to the DF. 

If we denote univariate LFs as cp^ (Li) and tjy^ {Lz), then the bivariate PDF (f>^ (Li, L2) is described by a differential copula 

c(u\ , U2) as 

(2) (L 1 ,L 2 ) = c[^ 1) (L 1 ),4 1) (L 2 )] . (32) 
For the FGM copula, the BLF leads from Equation i20\ 

<^\L^L 2 ;k) = {l + k [2*W(Li) - l] [29^ (L a ) - l] } <^\L{)^\L % ) . (33) 
The parameter k is proportional to the correlation coefficient p between log Li and log L2. For the Gaussian copula, the BLF is obtained as 

d>^(L 1 ,Lr,p) = ^=== exp |-i [*- lT (5T 1 - I) * ^] J ^(L^™ (I*) , (34) 

where 

= ($P(Li)) , *^ ($ 2 1) (i 2 ))] T (35) 

and S is again defined by Equation j25\ . 



3.2 The FIR-FUV BLF 



Here, to make our model BL F astrophysically reali stic, we construct the FUV-FIR BLF by the copula method. For the IR, we use the analytic 
form for the LF proposed by Saunders et alJ 1 1990h which is defined as 



L 



exp 



1 

2^2 



l0g ( 1 + I^ 



(36) 



We adopt the parameters estimated by Takeuchi et alJ j20 03bh which are obtained from the IRAS PSCz galaxies 1 Saunders et al. 2000). For 
the UV, we adopt the Schechter function 1 S^hechtelfl97^ ). 



4 1) (L) = (lnl0) <Mf^ 



exp 



L 
IT2 



(37) 



We use the parameters presented bv lWvder et al.U2005h for GALEXFUV (A„ ff = 1530 A): (a 2 , L t2 , ^.2) = (1.21, 1.81xl0 9 /i" 2 L , 1.35X 
10~ 2 h 3 Mpc -3 ). For simplicity, we neglect the A'-correction. We use the re-normalized version of Eqs. ( 1361 ) and {37} so that they can be 
regarded as PDFs, as mentioned above. 



3.3 Result 

We show the constructed BLFs from the FGM and Gaussian copulas in Figures [2] and [5] respectively. The FGM-based BLF cannot have a 
linear correlation coefficient larger than ~ 0.3 as explained above, while the Gaussian-based BLF may have a much higher linear correlation. 
We note that both copulas allow negative correlations, which are not discussed in this article. 

First, even if the linear correlation coefficients are the same, the detailed structures of the BLFs with the FGM and Gaussian copulas are 
different (see the case of p = 0.0-0.3). For the Gaussian-based BLFs, we see a decline at the faint end, while we do not have such structure 
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10* 10 s 10 10 10 

60/im luminosity L a0 [L Q ] 



Figure 3. The analytical BLF constructed with the Gaussian copula with model luminosity functions (LFs) of ultraviolet (UV) and infrared (IR)-selected 
galaxies.. The BLFs are again normalized so that integrating it over the whole ranges of L\ and L2 gives one. The linear correlation coefficient p varies from 
0.0 to 0.9 from top-left to bottom-right. The same as Fig. [2] the contours are logarithmic with an interval A log </>( 2 ' = 0.5 drawn from the peak probability. 



© 0000 RAS. MNRAS 000, 000-000 



Bivariate distribution function: application to the LF 9 



in the FGM-based BLFs (see the closed contours in Fig.O. This structure is introduced by the Gaussian copula, and from the physical point 
of view, it might not be strongly desired. The FGM-based BLF has a more ideal shape. 

Second, since the univariate LF shapes are different at FIR and FUV, the ridge of the BLF is not a straight line but clearly nonlinear. 
This feature is more clearly visible i n higher correlatio n cases in Figure[5] but always exists for the whole range of p. This trend is indeed 
found in t he Lfir-Lfuv diagram l IMartin et al . 2005). The underlying physics is that galaxies with high SFRs are more extinguished by 
dust (e.g. lBuat et al J2007al lbb. 

Observational applications including this topic will be presented elsewhere (Takeuchi et al. 2010, in preparation). 



4 DISCUSSION 

4.1 Flux selection effect in multiband surveys 

Since we have an explicit form of a BLF, we can discuss the flux selection effect formally. For simplicity, we consider the bivariate case 
(i.e. sample selected at two bands), but it will be straightforward to extend the formulation to a multiwavelength case (or more generally, 
selected using any physical properties). The flux selection is described in terms of luminosity as putting a lower bound L on a luminosity- 
luminosity (L1-L2) plane. The lower bound luminosity L llm is defined by the flux (density) detection limit S hrn as a function of redshift. In 
most surveys, a certain wavelength band is chosen as the primary selection band, like 6-band, ifs-band, 60 ^im-selected, etc. The schematic 
description of a survey is presented in Figure [4] 

If we select a sample of objects (in our case galaxies) at band 1, the objects with Li < Li (z) would not be included in the sample 
at a certain redshift z. Then, the detected sources should have Li > L\ m {z) and L2 > L 2 lm (z). Hence, on the L1-L2 plane, the 2-dim 
distribution of the detected objects is expressed as 



^-\dct 



(Li, La) = J' £^4> (2) (Li, L 2 ) 9 (Li im (z')) 6 (l?V)) & > (38) 
where Q, is a solid angle, and O is the Heaviside step function defined as 
fo L < a 

0(a)= . (39) 

[1 L ^ a 

The quantity S dct is proportional to the surface number density of objects detected in both bands on the L1—L2 plane. We start from a 
primary selection at band 1, then we would have objects detected at band 1 but not detected at band 2. In such a case we only have upper 
limits for these objects. The 2-dim distribution of the upper limits at band 2 is similarly formulated as 

^\ Li ,l 2 )^ £ £^^( Ll ,L 2 )e(L l r(z')) [i-e(L l r( z '))]dz' . m 

The superscript UL2 stands for "upper limit at band 2". In statistical terminology, the upper limit case, i.e. we know there is an object 
but we do only have the upper (or lower) limits of a certain quantity, is referred to as "censored". Though we can define the distribution 
E UL2 (Li, L2) by Eq. J40b . since the sample objects belonging to this category appear only as upper limits on the plot, a special statistical 
treatment, referred to as the survival analysis, is required to estimate E UL2 (Li, L2) from the data. Since we select objects at band 1, we 
do not have upper limits at band 1, because we do not know if there would be an object below the limit. This case is called "truncated" in 
statistics. 

If we select objects at band 2, we can formulate the 2-dim distribution of detected objects and upper limits exactly in the same way as 
the band 1 selected sample. For the objects detected at both bands, the 2-dim distribution is expressed by Eq. J38b . The objects detected at 
band 2 but not detected at band 1 is expressed as 



E UL1 ' 



(Lr,L 2 ) = f ^t^ (2> (Ll ' ia) I 1 - 9 ( L " m (^'))] 6 (^""V)) dz ' • C 41 ) 



If we can model L\ m (z) and L 2 lm (2) precisely including the isT-correction, evolutionary effect, etc., we can use the observed bivari- 
ate luminosity distribution to estimate the correlation coefficient, or more generally the dependence structure of two luminosities through 
Eqs. t38M40|. We can deal with these cases in a unified manner with techniques developed in survival analysis. We discuss this issue in a 
subsequent work (Takeuchi et al. 2010, in preparation). 



4.2 Other possible applications 

4.2.1 The star formation rate function 

The star formation rate (SFR) is one of the most fundamental quantities to investigate the formation and evolution of galaxies. The SFR is 
often estimated from the FUV flux of galaxies (or other related observables like Ha etc.) after "correcting" the dust extinction. However, 
some problems have been pointed out for this method. For instance, the relation between the UV slope ft (or equivalently, FUV— NUV 
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Figure 4. Schematic description of the selection effect in a multiband survey. Left: selection at band 1, center: selection at band 2, and right: bivariate selection 
at band 1 and 2. 



color) and the FIR-FUV flux ratio Lfir/Lfuv (often referred to as the IRX-/3 re lation) is frequently used to correct the extinction, but this 
relation is not always the same for various categories of star-forming galaxies (e.g. lBuat et al .120051 ; Boissi er et al .koollBoauien et alj20Qg|: 
Takeuchi et al koiOah . 

Instead, the total SFR obtained from the FUV and FIR luminos i ties would be a more reliable measure of the S F R since both are directl y 
observable values (e.g. llglesias-Paramo et al.l2004lBuat et al.l20oi ; |lglesias-Paramo et al.l2006l : lBuat et al.l2007al lbl jTakeuchi et al.l2010ah . 
Assuming a constant SFR over 10 8 yr, and Salpeter initial mass function jSalpeteilll955l . mass range: 0.1-100 Mq), we have the relation 
between the SFR and Lfuv 



log SFRfuv = log Z/fuv — 9-51 



(42) 



For the FIR, to transform the dust emission to the SFR, we assume that all the stellar light is absorbed by dust. Then, we obtain the following 
formula under the same assumption for both the SFR history and the IMF as those of the FUV, 

log SFRd ust = log Ltir - 9. 75 - log ( 1 - 77) . (43) 

Here, r\ is the fraction of the dust emission by old stars which is not related to the current SFR 1 Hirashita. Buat. & Inouell2003 ), and Ltir is 
the FIR luminosity integrated over A = 8-1000 /im. Thus, the total SFR is simply 

SFR tot = SFRfuv + SFR dust (44) 

(Igles ias-Paramo et alJbOQol) . Since the total SFR is basically estimated from the luminosities at FUV and FIR (note that SFRfuv oc Lfuv 
and SFRdust oc Ltir), the estimation of the PDF of the total SFR reduces to the estimation of the FIR-FUV BLF (e.g. lTakeuchi et alj2010bh . 
However, since the total SFR is the sum of two dependent variables, it is not straightforward to formulate the function unlike the cases we 
have seen above. This analysis will be discussed with a more specific methodology in our future work. 



4.2.2 The distribution of the specific star formation rate 



Another direct application is the distribution of the specific SFR (SSFR), SFR/M, where M* is the total stellar mass of a galaxy. The SSFR 
has gained much attention in the last decade, since the relation between M* and the SSFR of galaxies turns out to be a very important 
clue to understand the SF history of galaxie s: more massive galaxies have ceased their SF activity earlier in the cosmic time than less 
massiv e galaxies (downsizing in redshift: e.g. Cowie et"ai1ll99d : Boselli et alj|2001 ; Heavens et al.ll2004 ; Feulner et al. 2005 ; Noeske et al] 
l2007allbllPanter et alJl2007l;lDamen et alj|2009jbl among others). For a comprehensive summary of the downsizing, readers are encouraged 
to read Introduction of lFontanot et al.l ( 120091) . Despite of its importance, the treatment of multi wavelength data for this analysis is inevitably 
complicated and does not seem to be well understood, because we must deal with the data related to SFR and M, estimation. This might be, 
at least partially, the reason why the quantitative values of the M, -SSFR relation are different among different studies. 

As may easily guess after the above discussions, the M* -SSFR relation can be reduced to the relation between a luminosity at a certain 
mass-related band (often near IR bands) and a SF-related one (FUV, FIR, etc.) Then, we can model, for example, a Lk-Ltik bivariate 
luminosity function (Ly: A~-band luminosity jf| to examine the observed relation including all the selection effects. This is particularly useful 
for this topic, since iTakeuchi et all fcoiObl) found that the SFRF cannot be described by the Schechter function unlike the assumptions 
adopted in previous studies, but much more similar to the Saunders IR LF [Eq. 136H . The selection effect would be more complicated than 
in the case of the same Schechter marginals, but can be treated in the same way as discussed above. Thus, this may also be an interesting 
application of the copula-based BLF in the epoch of future large surveys. 



More precisely, L x -^tir-^fuv trivariate function might be appropriate for this issue. 
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5 SUMMARY AND CONCLUSIONS 

In this work, we introduced an analytic method to construct a bivariate distribution function (DF) with given marginal distributions and 
correlation coefficient, by making use of a convenient mathematical tool, called a copula. Using this mathematical tool, we presented an 
application to construct a bivariate LF of galaxies (BLF). Specifically, we focused on the FUV-FIR BLF, since these two luminosities are 
related to the star formation (SF) activity. Though both the FUV and FIR are related to the SF activity, the marginal univariate LFs have a 
very different functional form: former is well described by Schechter function whilst the latter has a much more extended power-law like 
luminous end. We constructed the FUV-FIR BLFs by the FGM and Gaussian copulas with different strength of correlation, and examined 
their statistical properties. Then, we discussed some further possible applications of the BLF: the problem of a multiband flux-limited sample 
selection, the construction of the SF rate (SFR) function, and the construction of the stellar mass of galaxies (A/*)-specific SFR (SFR/M,) 
relation. 

We summarize our conclusions as follows: 

(i) If the correlation of two variables is weak (Pearson's correlation coefficient \p\ < 1/3), the Farlie-Gumbel-Morgenstem (FGM) copula 
provides an intuitive and natural way for constructing such a bivariate DF. 

(ii) When the linear correlation is stronger, the FGM copula becomes inadequate, in which case a Gaussian copula should be preferred. 
The latter connects two marginals and is directly related to the linear correlation coefficient between two variables. 

(iii) Even if the linear correlation coefficient is the same, the structure of a BLF is different depending on the choice of a copula. Hence, a 
proper copula should be chosen for each case. 

(iv) The model FIR-FUV BLF was constructed. Since the functional shape of the LF at each w avelength is very different, the obtained 
BLF has a clear nonlinear structure. This feature was indeed found in actual observational data (e.g., Martin et alj2005h . 



(v) We formulated the problem of the multiwavelength selection effect by the BLF. This enables us to deal with datasets derived from 
surveys presenting complex selection functions. 

(vi) We discussed the estimation of the SFR function (SFRF) of galaxies. The copula-based BLF will be a convenient tool to extract 
detailed information from the observationally estimated SFRF because of its bivariate nature. 

(vii) The stellar mass-specific SFR relation was also discussed. This relation can be reduced to a BLF of luminosities at a mass-related 
band and a SF-related band. With an analytic BLF model constructed by a copula will provide us with a powerful tool to analyze the 
downsizing phenomenon with addressing the complicated selection effects. 

As the copula becomes better known to the astrophysical community and statisticians develop the copula functions, we envision many 
more interesting applications in the future. In a series of forthcoming papers, we will present more observationally-oriented applications of 
copulas. 
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APPENDIX A: AN EXTENSION OF THE FGM SYSTEM BY JOHNSON & KOTZ (1977) 



Although the FGM system of distributions provides us with convenient tool to construct the statistical model, its usefulness is restricted 
by the limitation of the corr elation s trength d escribed above. To overcome this dra wback, many attempt s have been mad e to extend the 
FGM distributions (see, e.g., Stuart et aHll994l : lKotz. Balakrishnan. & Johnsonll200ol) . Among them, Johnson & Kotz 1 1977) introduced the 
following iterated generalization of Equation il lb : 



(Al) 



G(x 1 ,x 2 ) = Y l K j [F 1 (x 1 )F 2 (x 2 )] [3/2]+1 {[1- F 1 (x 1 )][l- F 2 (x 2 )]} l(j+1)/2] , 

3=0 

where the symbol in the exponent [j/2] means the maximum natural number which does not exceed j/2. We set kq = 1. iHuang & Kotz 
( 1984) examined the dependence structure of Equation dAlt especially for the case of k = 2, and showed that the correlation can be stronger 
for these extension. In the case of the one-iteration family (k = 2), we have the DF as 



G(x 1 ,x 2 ) = F 1 (x 1 )F 2 (x 2 ) {1 + ki [1 - Fi(xi)] [1 - F 2 (x 2 )] + k 2 F 1 (x 1 )F 2 (x 2 ) [1 - Fi(xi)] [1 - F 2 (x 2 )]} . 
The corresponding PDF is 

g(xi,x 2 ) = fi(x 1 )f 2 (x2) {1 + Ki [2F!(xi) - 1] [2F 2 (x 2 ) - 1] + k 2 F 1 (xi)F 2 (x 2 ) [mfa) - 2] [3F 2 (x 2 ) - 2]} . 
Then, just the same as the case of the original FGM distribution (k = 1), we obtain the covariance 

Cov(aii, x 2 ) = j fj f(xi - X!)(x 2 - x 2 )f 1 (x 1 )f 2 (x 2 ) 

x {1 + ki [2F!(xi) - 1] [2F 2 (x 2 ) - 1] + K 2 F 1 {x 1 )F 2 {x 2 ) [3Fi(n) - 2] [3F 2 (x 2 ) - 2]}dx 1 dx 2 
= Kif x 1 f 1 (x 1 )[2F 1 (x 1 ) - l]dx 1 j x 2 f 2 (x 2 ) [2F 2 (x 2 ) - 1] dx 2 

+k 2 J x 1 f 1 (x 1 )F 1 (x 1 )[3F 1 (x 1 )~2]dx 1 J x 2 f 2 (x 2 )F 2 (x 2 )[3F 2 (x 2 )-2]dx 2 

= AlKl + A 2 K 2 . 

Huang & Kotz 1 1984 ) obtained the parameter space for «i and k 2 a^f] 



(A2) 



(A3) 



|«i| s: l , 

^1 + K 2 > — 1 



«2 ^ 



3-ft:i + x /3(l-Ki)(3+'^y 



They showed that, for a positive correlation, 

. Kl 3lK2 

p ^ • 

P 3 240 



(A4) 

(A5) 
(A6) 

(A7) 
(A8) 



Under these conditions, we have p ^ 0.5072, which is considerably better than 1/3. The BLF constructed with the first-order iterated FGM 
copula is expressed as 

^ (2) (L 1 ,L 2 ) = 0< 1) (L 1 )4 1) (i 2 ) 

x{l + «i [2^>[ 1) (L 1 )-l\ [2^ 1} (L 2 )-l] +k 2 ^ 1) (L 1 )4 1) (L 2 ) [3<E-< 1) (L 1 )-2] (L 2 ) - 2] } . (A9) 

The FIR-FUV BLF by Equation |A9| is shown in Figure lATI Clearly the dependence between the two luminosities is stronger than the original 
FGM-based BLF. However, now it is not intuitive nor straightforward to relate these two parameters of dependence ki and n 2 to the linear 
correlation coefficient. 



APPENDIX B: ESTIMATORS OF THE NONPARAMETRIC DEPENDENCE MEASURES p s AND r 



Here we present the estimators of the nonparametric measure of de pendence introdu ced in Section [2T2l More detailed derivation and properties 
of these nonparametric measures of dependence are found in e.g., Hettmanspergej Jl984ri and lHollander & Wolf j Jl999ri . 

Let {Ri}i=i,...,n and { Si}i=i, be the ranks of {xii}j=i,..., n and {a; 2 i}i=i,...,n, respectively. If we denote the estimator of ps for a 
sample as rs, 

EIU m (Si - H 1 ) 



r s = 



^eiu (Ri - m Er=i (Si - s±i) 



(Bl) 



3 Not e that there is a typo for Equation (A6) in Huang & Kotz ( 1984). There are also typos in Equations jA6\ and )A7t in Kotz, Balakrishnan, & Johnson 
<200d). 
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Figure Al. The BLF constructed with a first-order iterated FGM copula. The linear correlation coefficient p is 0.45. 



This is also expressed as 

n(n z — 1) 

This form is an exact sample counterpart of Equation {8]l. If we define di = Si — -Ri, Equation (B2\ reduces to the following simper form 

6E" i dj 

rS = 1 ~ M ^ • (B3) 

The variance of rs in the large-sample limit is given by 

Var[r s ] = — l — . (B4) 
n — 1 

The most basic form of Kendall's r for a sample has been already shown in Section[2j2] [Equation J5J]. It is also expressed as 

2(n c — n d ) 
n(n — 1) 

= 1 - , 4 7r n d (B5) 
n(n — 1) 

since n c + rid = n(n — 1) /2. The variance of t is given by 

Var M = f^±§ (B6) 

971(71 — 1) 



(see IValz & Mcleodll 199(1 for a very concise derivation). 
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